home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / a2_7.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.7 KB  |  168 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 2.7 (Steffensen's Acceleration).
  9. % Section    2.5,    Aitken's Process & Steffensen's & Muller's Methods, Page 96
  10. echo on; clc; format long; hold off; clear
  11. % This program implements the Steffensen method.
  12.  
  13. %    Define and store the functions f(x) and f'(x)
  14.  
  15. %    in the M-files  f.m  and  df.m, respectively.
  16. % function y = f(x)
  17. % y = x.^3 - 3.*x + 2;
  18.  
  19. % function y1 = df(x)
  20. % y1 = 3*x.^2 - 3;
  21.  
  22. delete f.m
  23. diary f.m; disp('function y = f(x)');...
  24.            disp('y = x.^3 - 3.*x + 2;');...
  25. diary off;
  26.  
  27. delete df.m
  28. diary df.m; disp('function y1 = df(x)');...
  29.             disp('y1 = 3*x.^2 - 3;');...
  30. diary off;
  31.  
  32. % Remark. f.m, df.m and steff.m are used for Algorithm 2.7
  33. f(0); df(0); % Test for files f.m, df.m
  34. pause    % Press any key to see the graph y = f(x).
  35.  
  36. clc; clg;
  37. %    Plot f(x) over the interval [a,b].
  38.  
  39. a = -3.0;
  40. b = 3.0;
  41. c = -10;
  42. d = 10;
  43. h = (b-a)/150;
  44. X = a:h:b;
  45. Y = f(X);
  46. axis([a b c d]);...
  47. plot(X,Y,'-g');...
  48. hold on;...
  49. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  50. xlabel('x');...
  51. ylabel('y');...
  52. title('Graph of y = f(x).');...
  53. grid;...
  54. axis;...
  55. hold off;...
  56. shg; pause    % Press any key to perform Steffensen iteration.
  57.  
  58. clc;
  59. %    Place the starting value in  p0
  60.  
  61. %    Place the abscissa tolerance in  delta
  62.  
  63. %    Place the ordinate tolerance in  epsilon
  64.  
  65. %    Place the number of iterations in  max1
  66.  
  67. p0  = -2.4;
  68. delta = 1e-12;
  69. epsilon = 1e-12;
  70. max1 = 5;
  71.  
  72. [p,yp,err,Q] = steff('f','df',p0,delta,epsilon,max1);
  73.  
  74. pause    % Press any key for the list of iterations.
  75.  
  76. clc; clg;
  77. a = -2.42;
  78. b = -1.97;
  79. c = -5.0;
  80. d = 0.5;
  81. h = (b-a)/150;
  82. X = a:h:b;
  83. Y = f(X);
  84. max1 = length(Q);
  85. n0 = min(6,max1);
  86. X0 = Q(1:n0);
  87. Z0 = zeros(1,n0);
  88. axis([a b c d]);...
  89. plot(X,Y,'-g',X0,Z0,'or');...
  90. hold on;...
  91. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  92. xlabel('x');...
  93. ylabel('y');...
  94. title('Graphical analysis for Steffensen`s method.');...
  95. grid;...
  96. axis;...
  97. hold off;...
  98. shg; pause    % Press any key to continue.
  99.  
  100. J = 1:max1;
  101. Yq = f(Q);
  102. points = [J;Q;Yq];
  103. Mx1 = 'Iterations for Steffensen`s method.';
  104. Mx2 = ['     k                  p(k)               f(p(k))'];
  105. Mx3 = 'The solution is:';
  106. Mx4 = 'The error estimate for p is  ± ';
  107. clc,echo off,diary output,
  108. disp(''), disp(Mx1),disp(''), disp(Mx2), disp(points'),...
  109. disp('Iteration converged quadratically to the root.'),...
  110. disp(''),disp(Mx3),disp(''),disp('p = '),...
  111. disp(p),disp(''),disp('f(p) = '),disp(yp),...
  112. disp([Mx4,num2str(err)]),diary off,echo on
  113.  
  114. pause    % Press any key to perform Steffensen iteration.
  115.  
  116. clc;
  117. %    Place the starting value in  p0
  118.  
  119. %    Place the abscissa tolerance in  delta
  120.  
  121. %    Place the ordinate tolerance in  epsilon
  122.  
  123. %    Place the number of iterations in  Max
  124.  
  125. p0  = 1.2;
  126. delta = 1e-12;
  127. epsilon = 1e-12;
  128. max1 = 5;
  129.  
  130. [p,yp,err,Q] = steff('f','df',p0,delta,epsilon,max1);
  131.  
  132. pause    % Press any key for the list of iterations.
  133.  
  134. clc; clg;
  135. a = 0.975;
  136. b = 1.225;
  137. c = -0.02;
  138. d = 0.14;
  139. h = (b-a)/150;
  140. X = a:h:b;
  141. Y = f(X);
  142. max1 = length(Q);
  143. n0 = min(6,max1);
  144. X0 = Q(1:n0);
  145. Z0 = zeros(1,n0);
  146. axis([a b c d]);...
  147. plot(X,Y,'-g',X0,Z0,'or');...
  148. hold on;...
  149. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  150. xlabel('x');...
  151. ylabel('y');...
  152. title('Graphical analysis for Steffensen`s method.');...
  153. grid;...
  154. axis;...
  155. hold off;...
  156. shg; pause    % Press any key to continue.
  157.  
  158. J = 1:max1;
  159. Yq = f(Q);
  160. points = [J;Q;Yq];
  161. clc,echo off,diary output,
  162. disp(''), disp(Mx1),disp(''), disp(Mx2), disp(points'),...
  163. disp('Iteration converged quadratically to the root.'),...
  164. disp(''),disp(Mx3),disp(''),disp('p = '),...
  165. disp(p),disp(''),disp('f(p) = '),disp(yp),...
  166. disp([Mx4,num2str(err)]),diary off,echo on
  167.  
  168.